Theory of Branching and Annihilating Random Walks 
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A systematic theory for the diffusion-limited reaction processes A + A — > and A —* (m + I) A 
is developed. Fluctuations are taken into account via the field-theoretic dynamical renormalization 
group. For m even the mean field rate equation, which predicts only an active phase, remains 
qualitatively correct near d c = 2 dimensions; but below d' c ~ 4/3 a nontrivial transition to an 
inactive phase governed by power law behavior appears. For m odd there is a dynamic phase 
transition for any d < 2 which is described by the directed percolation universality class. 



PACS numbers: 64.60.Ak, 05.40+j, 64.60. Ht 

Nonequilibrium models with an extensive number of 
degrees of freedom whose dynamics violates detailed bal- 
ance occur in studies of many biological, chemical and 
physical systems. Like equilibrium systems, their station- 
ary states may exhibit phase transitions which in many 
cases appear to fall into distinct classes characterized by 
universal quantities such as critical exponents. One of 
the most common such classes is that exemplified by di- 
rected percolation (DP) This represents a transition 
from a nontrivial 'active' steady state to an absorbing 'in- 
active' state with no fluctuations. Many nonequilibrium 
phase transitions appear to belong to this universality 
class, e.g., the contact process [0, the dimer poisoning 
problem in the ZGB model ||, and auto-catalytic re- 
action models 0. The universal properties of the DP 
transition are theoretically well understood in the con- 
text of a renormalization group (RG) analysis based on 
an expansion around mean held theory below the upper 
critical dimension d c = 4 Q . 

More recently a class of models has been studied which, 
in certain cases, appear as exceptions to the general rule 
that such transitions should fall into the DP universal- 
ity class. These include a probabilistic cellular automa- 
ton model ||, certain kinetic Ising models and an 
interacting monomer-dimer model H. In one dimen- 
sion the dynamics of these is equivalent to a class of 
models called branching and annihilating random walks 
(BARWs) [|1q| — P~2f] 7 which also have a natural generaliza- 
tion to higher dimensions. In the language of reaction- 
diffusion systems, BARWs describe the stochastic dy- 
namics of a single species of particles A undergoing three 
basic processes: diffusion, often modeled by a random 
walk on a lattice and characterized by a diffusion coef- 
ficient D; an annihilation reaction A + A — > when 
particles are close (or on the same site), at rate A; and 
a branching process A — > (to + 1)^4 (where to is a pos- 
itive integer), at rate o~ m . The above-mentioned one- 
dimensional models all correspond to the case m — 2. For 
the kinetic Ising model, the particles A are to be iden- 
tified with the domain walls, and the transition to the 
inactive state corresponds to the ordering of the Ising 



spins In general, this new universality class has 

been observed in d = 1 for even values of m, when the 
number of particles is locally conserved modulo 2. When 
to is odd, the DP values of the exponents appear to be re- 
alized. (It should be remarked that several of the models 
which have been studied do not contain three indepen- 
dent parameters corresponding to D, A, and o~ m so that it 
may occur that the actual transition is inaccessible. This 
appears to be so for the simplest lattice BARW model 
with to = 2, which is always in the inactive phase [fl0||.) 

Besides the appearance of a new universality class, an- 
other issue which clearly requires theoretical explanation 
is the occurrence of a transition at a finite value of o~ m . 
For the mean held rate equation for the average density 



h(t) = -2\n(t) 2 + 



ma, 



n(t) 



(1) 



predicts a non-zero steady state density mo~ m /2\, so that 
this state should be active for all o~ m > 0, in contrast to 
what is in fact observed for d = 1. It is, however, known 
that fluctuation effects render the mean field description 
of the pure annihilation problem (<r m = 0) qualitatively 
incorrect for d < 2 Jf3|. Therefore, a detailed theory has 
to demonstrate how these are also responsible for moving 
the critical value of o~ m away from zero. 

In this Letter, we describe the first systematic theory 
of these phenomena (details will be presented elsewhere 
[|l4||). It is based on the field-theoretic RG analysis which 
has proven successful for the DP problem pj: however, 
we correct an important error which was made in an ear- 
lier investigation along the same lines jq] . A summary of 
our main results follows. 

(a) Even to. For d > 2 the mean field description (|I|) 
is qualitatively correct in that the transition occurs at 
o~ m = 0. The density in the active phase vanishes as 
n oc o~ m , with calculable logarithmic corrections in two 
dimensions. As d is lowered below 2, the transition first 
continues to occur at o~ m = 0, with modified critical ex- 
ponents, until a second critical dimension d' c ks 4/3 is 
reached. Below this, and in particular for d = 1, there 
appears a nontrivial transition at o~ c > from the active 
phase to an inactive phase in which the density decays 
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asymptotically as n oc t~ d / 2 . Because of the existence of 
two critical dimensionalities, this new universality class 
apparently has no simple mean field limit, close to which 
the fluctuations can be controlled. We are therefore un- 
able to generate a systematic e expansion for the critical 
exponents. The truncated loop expansion in fixed dimen- 
sion seems to provide at least a qualitative description of 
the transition; however, as there exists no small expan- 
sion parameter, the actual values for the critical expo- 
nents to one-loop order are rather inaccurate. The RG 
analysis also shows that higher values of m inevitably 
generate an effective m = 2 reaction under renormaliza- 
tion, and this is always the most relevant term. Therefore 
all such processes with even m fall into the same univer- 
sality class. We have also considered an N species gen- 
eralization of the rn = 2 model. The RG analysis shows 
that for N > 1 the critical behavior is that of the N — » oo 
limit, which is exactly solvable but always in the active 
phase for ai > 0, with critical exponents being described 
by yet another universality class. 

(b) Odd m. The case m = 1 here is typical. Under 
renormalization, a spontaneous decay process A — » 
is generated for arbitrarily small branching rates, by the 
combined reactions A — > 2 A, 2 A — > 0, so that the density 
decays exponentially, as in the inactive phase of DP. We 
find this to occur for d < 2, when the pure annihilation 
process is relevant. For d > 2, however, the situation is 
more subtle. The RG analysis in this case predicts that 
there is in fact a nontrivial transition even at d c = 2, 
with ct c ~ De~ i7tD / x , while it is absent for d > 2. Anal- 
ysis of the effective theory for d < 2 then shows that the 
subsequent transition at larger values of o\ is in the DP 
universality class, as is observed in simulations jl(|[ll| • 

The field-theoretic analysis of these problems be- 
gins from the 'second-quantized' approach to classical 
stochastic particle systems which is well known and 
has been described in detail elsewhere EJl^l . Anni- 
hilation and creation operators a* and a], satisfying 
the usual boson commutation relations, are introduced 
at each site i of the lattice, and the time-dependent 

state vector = J2{n i }P({ n i}' t )Hi a i is con " 

structed from the probabilities p({ni};t). The (clas- 
sical) master equation satisfied by these may then be 
recast as a Schrodinger— like equation with a time evo- 
lution operator which, in this example, has the form 
H = H d + H a + H b + H T , where 

(2) 
(3) 
(4) 
(5) 



to) 

i 

Ef 4-m+l + 
I a\ di - a\a 

i 

H T = -rJ2( a f- 1 



The last term corresponds to a constant creation of 
pairs of particles, simulating the effects of finite temper- 
ature in the kinetic Ising model ||. Finally, utilizing 
the coherent-state path integral formalism, the 'quan- 
tum many particle' hamiltonian H can be cast into a 
field theory which describes BARW processes including 
fluctuation effects. Note that no additional assumptions, 
specifically regarding the form of the noise correlations in 
an extension of Eq. ([!]) to an effective Langevin equation, 
had to be invoked in this derivation. 



When m is even, there is a formal symmetry of H 
under changing the signs of all the <n and a\ simulta- 
neously: this corresponds to the conservation of parti- 
cle number modulo 2. However, in this formalism, ex- 
pectation values of operators such as the local density 
= a\aj, for example, are given by matrix elements 



71 j — ujuj 

of the form 



Ht 



^(0)), and in order to use 



time-dependent perturbation theory and Wick's theorem 

it is conventional to commute the factor e^-^ a * through. 
This is equivalent to applying the formal shift at — > 1+a] 
in H . Yet this obscures the above symmetry, and if, in ac- 
cordance with the usual naive power counting arguments 
near the upper critical dimension, higher order quartic 
terms in H are then ignored, it becomes completely lost. 
This led the authors of Ref. || to the erroneous conclu- 
sion that, near d = 4, the transition should be in the DP 
universality class irrespective of the parity of m. 

However, it is imperative in any RG analysis to pre- 
serve all known symmetries of the system. In the present 
case, this may be done by observing that the RG equa- 
tions themselves (as opposed to the calculations of ob- 
servables such as the density) should be independent of 
which basis is used, and it is therefore possible, and, in- 
deed, necessary, to perform the computations in the rep- 
resentation of the model in which the symmetry is mani- 
fest. The methods for doing this are standard, and will be 
described in detail elsewhere 14 1. The case a m = r = 0, 
corresponding to a pure annihilation process, has already 
been analyzed in [^3|. The RG equation for the flow of 
the dimcnsionlcss coupling I = CdX/Dn^, where n is a 
normalization wave number, Cd = T(2 — d/2)/2 d_1 7r d / 2 
a geometric factor, and e = 2 — d, under a rescaling factor 
e, is given by dt/dl — e£ — £ 2 , which is exact at one loop. 
For d < 2 the late time behavior is controlled by the 
nontrivial fixed point at £* — e, leading to an asymptotic 
particle density decay according to n(t) oc t~ d f 2 . 

The first question to be addressed is whether the 
branching rate a m is relevant at the pure annihilation 
fixed point, i.e., whether its RG eigenvalue y a is posi- 
tive. If so, the late time behavior must differ from that 
of the pure annihilation process, indicating that the ac- 
tive phase is reached immediately. For d > 2 we find 
y a — 2 from simple power counting, so indeed o~ m is rele- 
vant. The density in the active phase vanishes according 



2 



to the mean field result n oc a m . For d < 2, to one-loop 
order, y a = 2 — {m(m+l)/2]£+0(£ 2 ), so that a m remains 
relevant at the annihilation fixed point £* = e just below 
d — 2, with the lower values of m the being most relevant. 
Since these lower allowed values of to inevitably become 
generated whenever the annihilation rate is nonzero, we 
conclude that the cases with to even will always fall into 
the universality class of to = 2, while to odd will generate 
m = 1 and to = — 1. The latter is always relevant, and, 
as we shall see below, is responsible for the crossover to 
the DP universality class. For the time being we there- 
fore restrict our attention to the case of even to. For 
d = 2 the marginality of I is responsible for logarithmic 
corrections to mean field theory, which for m = 2 take 
the form n oc cr/[ln(l/cr)] 2 . 

The above result for y a is valid only close to d = 2. For- 
tunately it is possible to compute it exactly in one dimen- 
sion, at the pure annihilation fixed point. The latter cor- 
responds to the limit of infinite bare coupling A |ll| . The 
multiparticle states then effectively propagate as hard- 
core bosons in between the annihilation and branching 
processes, and so, in one dimension, behave like free 
fcrmions. On the lattice, this limit only makes sense if we 
define the branching process as placing the to offspring on 
different but neighboring sites. The branching contribu- 
tion to H, in terms of these fermionic operators c, and cj, 

thus acquires the form H b = a m £\ n^f^ TO / 2 c i+j°i- Tlic 
continuum limit of this expression, found by performing 
a Taylor expansion in powers of the lattice spacing a , 
will be different from the bosonic case because the an- 
ticommuting nature of the c\ allows each derivative to 
appear only once. The lowest order term has the form 
a™ (m+1)/2 ct(<9ct)(d 2 ct) . . . (<9 m ct)c, with the result that 
the effective expansion parameter a m = a ™( m + 1 )/ 2 '• (Tm 
has a modified scaling dimension. This leads to the result 
(which may be confirmed by other less formal methods) 
that y a = 2 — to (to + l)/2 exactly in d = 1. Thus, for 
reasons we do not understand, the 0(e) result appears 
to be exact in d = 1, and y a changes sign at a value 
of d = d' c w 4/3 for m — 2, if the higher order terms 
continue to be small. In d = 1, y CT < for all the even 
values of to. This establishes the result that the late time 
behavior for small values of a m is controlled by the anni- 
hilation fixed point, so that n(t) oc t~ x l 2 . In the inactive 
phase, the system is composed of a set of highly anticor- 
related bunches of odd numbers of particles, the spatial 
distribution of which, upon coarse graining, looks like 
that of single particles in the pure annihilation process. 

Clearly, the above scenario cannot be obtained in any 
finite order of an expansion near d c — 2. We have there- 
fore performed a truncated loop expansion at fixed di- 
mension, retaining the full dependence on a, which ap- 
pears both as a vertex and as a mass term. To one loop 
order, the RG flow equations for the renormalized reac- 
tion rates I = Cd\/DK 2 ~ d and s = a/Dn 2 read (m = 2) 



d£/dl = £[2-d- 1/(1 + sf- d l 2 ] , (6) 
ds/dl = s[2 - 3£/(l + s) 2 - d ' 2 ] . (7) 

For s — > 0, the annihilation fixed point £* = 2 — d of the 
inactive phase is recovered, while for s — > oo the loop con- 
tributions to the anomalous dimensions vanish, and the 
flow approaches a Gaussian fixed point describing the ac- 
tive phase. For large s, the effective coupling in Eqs. (||), 
(0) becomes g = £/s 2 ~ d ' 2 (see Sec. Ill of Ref. @), whose 
flow is given by dg/dl = 2g- [(10-3c?)/2] 5 2 = -(3(g). In 
addition to the stable Gaussian fixed point at g = there 
is a nontrivial unstable one at g* = 4/(10 — 3c?) describing 
the phase transition. At this order there is neither field 
nor diffusion constant renormalization, giving a dynamic 
exponent z ~ 2. However, because the mean field density 
n ~ u/A and the spatial correlation length £j_ ~ cr -1 / 2 
depend not just on g but also on the dangerous irrele- 
vant variable s _1 , the critical exponents describing the 
approach to the critical point in the active phase depend 
not only on y E = (3'(g*), describing the distance from the 
critical point s = (g* — g) /g*, but also on y\ w 2 — d — g* 
and y a « 2 — 3g*. As a result we find that n ~ e 13 , with 
P = (d + yx- y*)/Ve « 4/(10 - 3d), and £ ± - e"^, 
with v± = (1 — y a )/ye ~ 3/(10 — 3c?). The truncated 
one-loop approximation thus seems to provide a qual- 
itatively correct picture of the transition, although the 
actual numerical values of these exponents in one dimen- 
sion are rather poor as compared to simulation results 
this is not very surprising, however, as there is no 
small expansion parameter present here. In addition, wc 
cannot really access those exponents that describe the 
behavior at the critical point, as the density might de- 
pend nonanalytically on a there; this also precludes a 
sound derivation of scaling relations || connecting these 
with the above exponents describing the active phase. 

A better result is obtained for the exponent v T = l/y T 
describing the divergence of the correlation length as the 
pair creation rate r — > at the critical point e = 0. The 
one-loop flow equation for r reads 

dr/dl = r[d + 2- £/(! + sf- d l 2 \ , (8) 

and hence in the inactive phase, or at the critical point 
a = ior d > d' c , one has v T = l/2d, while at the non- 
trivial phase transition for d < d' c the result is v r w 
(10 - 3d)/(16 + Ad - 3o! 2 ). In one dimension, v r w 7/17, 
which is in fair agreement with simulations S. 

We have also investigated an N species generalization 
of the to = 2 problem, defined by the processes 2A a — * 0, 
at rate A/A, A a -> 3A a , at rate a, and A a —> A a + 2A 13 , 
(3 ^ a, at rate cr'/(N — 1). To one-loop order at the 
annihilation fixed point the RG eigenvalue of the addi- 
tional branching process becomes y a > = 2 — £, which is 
therefore more relevant than the original reaction with 
rate a. We have chosen the above N component version, 
because for N — > oo, the ensuing theory (with a = 0) can 
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be solved exactly; physically this limit corresponds to the 
situation where each particle may annihilate only with its 
sibling. The resulting critical point remains at a c = for 
all d, and its universality class, distinct from the previ- 
ously discussed ones, is characterized by the mean-field 
exponents z = 2, /3 = 1, and as a consequence of the now 
exact result y a i — 2 — £, we find v±_ = 1 = 1 /d, using 
t = 2 - d for d < 2. 

We now return to the case of odd m. As argued above, 
fluctuations generate a spontaneous single-particle de- 
cay process, and the effective interactions at a given site 
acquire the form 



(a) -l) a -a (a f - l) a) a + A (a f2 - l) a 2 



(9) 



When the single-particle decay rate \i ^ 0, it is con- 
venient to remove the linear term in a by the shift 
a' — > 1 + <v mentioned earlier. This results in the in- 
teraction hamiltonian 

(/i — a) a) a — a (r 



rrmt 

■"cff 



a + 2Aa t a 2 + Aa t a 2 . (10) 



If we now neglect the quartic term (justifiable in this 
case since there is no 'parity' symmetry that must be 
preserved), we find precisely the interaction hamiltonian 
used to characterize DP ||. The transition occurs when 
the renormalized version of the mass term \ir — an van- 
ishes. The question of whether this actually happens for 
allowed values of the bare parameters a and A depends 
on the way these are renormalized, and this may be stud- 
ied close to d — 2. It is simpler to work in the unshifted 
version (0), where it is straightforward to identify the 
most singular ('bubble') diagrams in powers of e _1 at a 
given order in A. The mass in the shifted DP hamiltonian 
( |l0| ) then becomes fin — an — a(Jd — l)/(-f<2 + 1); where 
I d = (C d X/De)[( f i R + a R )/D}-^ 2 = (C d \/ De){a / D)^/ 2 
because [Ir + a R = a in this approximation. For suffi- 
ciently small <7, this is positive, indicating that the sys- 
tem is in the inactive phase with an exponential decay 
of the density. The transition to the active phase occurs 
at a c = D{Cd\/ 'De) 2 / e . Although this result is accurate 
only to leading order in e, the general feature of a tran- 
sition at a finite value of a in the DP universality class 
should persist to d = 1. For d = 2 the transition is seen 
to continue to occur at a finite value a c ~ Der^ D l x , 
as A — > 0. However, for larger d, the annihilation rate 
A which drives the generation of the process A — > 0, es- 
sential for the DP inactive state, becomes irrelevant, and 
one may use the same set of diagrams to argue that there 
is now no transition, at least for small A. 

The same result can be obtained in the framework of an 
RG calculation similar to that invoked for m even. This 
method yields for m — 3 the same qualitative picture as 
for m = 1, but the transition moves closer to the mean 
field critical point, with a c ~ De~ 5 - 68nD / x Q, and we 
expect this tendency to hold for larger odd values of m, 
in accord with numerical simulations [IlOllLlf . 



To summarize, we have provided the first analytic the- 
ory of branching and annihilating random walks which 
explains most of their observed behavior. We have shown 
that the fluctuations responsible for the failure of mean 
field theory in the pure annihilation process for d < 2 
are also responsible for shifting the critical value of the 
branching rate away from zero. For m odd this occurs 
for all d < 2, with the subsequent transition being in 
the DP universality class, while for m even this effect is 
postponed to lower dimensions d < d' c « 4/3. Our theory 
correctly takes account of the symmetry in this case, but 
is so far unable to yield accurate estimates of the critical 
exponents in d = 1. However, a truncated loop expan- 
sion appears to provide at least a qualitative picture of 
the transition. It would of course be desirable to find 
some other controlled approximation scheme in which to 
approach this problem. Our investigation of an N species 
generalization of the m = 2 reaction failed to provide us 
with additional insight in the single species case, but in- 
stead uncovered yet another new universality class for 
N > 1, with a c = and governed by the exponents 
of the exactly solvable N — > oo limit. This underlines 
the importance of fluctuations and correlation effects in 
reaction-diffusion systems at low dimensions, which may 
lead to remarkably rich nonequilibrium phase diagrams. 
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